/*
 * RELIC is an Efficient LIbrary for Cryptography
 * Copyright (C) 2007-2013 RELIC Authors
 *
 * This file is part of RELIC. RELIC is legal property of its developers,
 * whose names are not listed here. Please refer to the COPYRIGHT file
 * for contact information.
 *
 * RELIC is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * RELIC is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with RELIC. If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * @file
 *
 * Implementation of multiple precision operations in polynomial rings.
 *
 * @version $Id: relic_bn_pol.c 1542 2013-08-28 18:31:47Z dfaranha $
 * @ingroup bn
 */
 

#include "relic_core.h"

/*============================================================================*/
/* Private definitions                                                        */
/*============================================================================*/

/**
 * Size of precomputation table.
 */
#define TABLE_SIZE			64

/*============================================================================*/
/* Public definitions                                                         */
/*============================================================================*/

void pol_init(pol_t a, dig_t deg) {
	a->degree = deg;
	a->coef = (bn_t *) calloc(deg+1, sizeof(bn_t));
}

void pol_zero(pol_t a) {
	int i;
	for (i=0; i<=a->degree; i++) {
		bn_zero(a->coef[i]);
	}
}

void pol_null(pol_t a) {
	int i;
	for (i=0; i<=a->degree; i++) {
		bn_null(a->coef[i]);
	}
}

void pol_new(pol_t a) {
	int i;
	
	TRY {
		for (i=0; i<=a->degree; i++) {
			bn_new(a->coef[i]);
		}
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
}

void pol_free(pol_t a) {
	int i;
	
	for (i=0; i<=a->degree; i++) {
		bn_free(a->coef[i]);
	}
	free(a->coef);
	free(a);
}

void pol_copy(pol_t c, pol_t a) {
	int i;
	
	pol_init(c, a->degree);
	for (i=0; i<=a->degree; i++) {
		bn_copy(c->coef[i], c->coef[i]);
	}
}

int pol_degree(pol_t a) {
	return a->degree;
}

void pol_coef(bn_t c, pol_t a, dig_t i) {
	bn_copy(c, a->coef[i]);
}

void pol_set_coef(pol_t a, dig_t i, bn_t value) {
	if (i<=a->degree) {
		bn_copy(a->coef[i], value);
	}
}

void pol_trim(pol_t a) {
	int i;
	
	for (i=a->degree; i>=0; i--) {
		if (!bn_is_zero(a->coef[i])) {
			break;
		}
	}
	
	if (a->degree != i) {
		a->coef = (bn_t *) realloc(a->coef, sizeof(bn_t)*(i+1));
		a->degree = i;
	}
}


void pol_print(pol_t a) {
	int i;
	for (i=0; i<=a->degree; i++) {
		printf("a%d: ", i);
		bn_print(a->coef[i]);
	}
	
}

int pol_cmp(pol_t a, pol_t b) {
	int i, cmp;
	
	if (a->degree > b->degree) {
		return CMP_GT;
	}
	else if (a->degree < b->degree) {
		return CMP_LT;
	}
	for (i=a->degree; i>=0; i--) {
		cmp = bn_cmp(a->coef[i], b->coef[i]);
		if (cmp != CMP_EQ) {
			return cmp;
		}
	}
	return CMP_EQ;
}

void pol_add(pol_t c, pol_t a, pol_t b) {
	int i, max=2;
	
	//max = (int) fmax((double)a->degree, (double)b->degree);
	pol_init(c, max);
	
	for (i=0; i<=max; i++) {
		if (i > a->degree) {
			bn_copy(c->coef[i], b->coef[i]);
		}
		else if (i > b->degree) {
			bn_copy(c->coef[i], a->coef[i]);
		}
		else {
			bn_add(c->coef[i], a->coef[i], b->coef[i]);
		}
	}
}


void pol_sub(pol_t c, pol_t a, pol_t b) {
	int i, max=2;
	
	//max = (int) fmax((double)a->degree, (double)b->degree);
	pol_init(c, max);
	
	for (i=0; i<=max; i++) {
		if (i > a->degree) {
			bn_copy(c->coef[i], b->coef[i]);
		}
		else if (i > b->degree) {
			bn_copy(c->coef[i], a->coef[i]);
		}
		else {
			bn_sub(c->coef[i], a->coef[i], b->coef[i]);
		}
	}
}










/**---------------------------------------------------------------*/








void bn_pol_null(bn_t *c) {
	bn_null(c[0]);
	bn_null(c[1]);
}

void bn_pol_new(bn_t *c) {
	TRY {
		bn_new(c[0]);
		bn_new(c[1]);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
}

void bn_pol_free(bn_t *c) {
	bn_free(c[0]);
	bn_free(c[1]);
}

void bn_pol_copy(bn_t *c, bn_t *a) {
	bn_copy(c[0], a[0]);
	bn_copy(c[1], a[1]);
}

void bn_pol_sqr_u(bn_t *c, bn_t *a, bn_t *m, bn_t n, bn_t u) {
	bn_t t0, t1, a0a1, r;
	
	bn_null(t0);
	bn_null(t1);
	bn_null(a0a1);
	bn_null(r);
	
	TRY {
		bn_new(a0a1);
		bn_new(t0);
		bn_new(t1);
		bn_new(r);
		
		bn_mul(a0a1, a[0], a[1]);
		bn_mod(a0a1, a0a1, n, u);
		
		bn_add(t0, a[0], a[1]);
		if (bn_cmp(t0, n) != CMP_LT) {
			bn_sub(t0, t0, n);
		}
		
		bn_pol_mul_swt(t1, a[1], m[0], n);
		//bn_mul(t1, a[1], m[0]);
		//bn_mod(t1, t1, n, u);
		
		bn_add(t1, t1, a[0]);
		if (bn_cmp(t1, n) != CMP_LT) {
			bn_sub(t1, t1, n);
		}
		
		bn_mul(t0, t0, t1);
		bn_mod(t0, t0, n, u);
		
		bn_pol_mul_swt(t1, a0a1, m[0], n);
		//bn_mul(t1, a0a1, m[0]);
		//bn_mod(t1, t1, n, u);
		
		bn_sub(t0, t0, t1);
		if (t0->sign == BN_NEG) {
			bn_add(t0, t0, n);
		}
		bn_sub(c[0], t0, a0a1);
		if (c[0]->sign == BN_NEG) {
			bn_add(c[0], c[0], n);
		}
		
		bn_lsh(c[1], a0a1, 1);
		if (bn_cmp(c[1], n) != CMP_LT) bn_sub(c[1], c[1], n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(a0a1);
		bn_free(r);
	}
}


void bn_pol_sqr_u_old(bn_t *c, bn_t *a, bn_t *m, bn_t n, bn_t u) {
	bn_t t0, t1, r;
	
	bn_null(t0);
	bn_null(t1);
	bn_null(r);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(r);
		
		bn_sqr(r, a[1]);
		bn_mod(r, r, n, u);
		
		bn_mul(r, r, m[0]);
		bn_mod(r, r, n, u);
		
		bn_mul(t1, a[0], a[1]);
		bn_mod(t1, t1, n, u);
		bn_lsh(c[1], t1, 1);
		if (bn_cmp(c[1], n) != CMP_LT) {
			bn_sub(c[1], c[1], n);
		}
		
		bn_sqr(t0, a[0]);
		bn_mod(t0, t0, n, u);
		
		bn_add(c[0], t0, r);
		if (bn_cmp(c[0], n) != CMP_LT) {
			bn_sub(c[0], c[0], n);
		}

	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(r);
	}
}

void bn_pol_sqr_basic(bn_t *c, bn_t *a, bn_t *m, bn_t n) {
	bn_t t0, t1, r;
	
	bn_null(t0);
	bn_null(t1);
	bn_null(r);
	
	TRY {
		bn_new(r);
		bn_new(t0);
		bn_new(t1);

		bn_sqr(r, a[1]);
		bn_mod(r, r, n);
		
		//bn_mul(r, r, m[0]);
		bn_pol_mul_swt(r, r, m[0], n);
		
		bn_sqr(t0, a[0]);
		bn_add(t0, t0, r);
		bn_mul(t1, a[0], a[1]);
		bn_lsh(t1, t1, 1);
		
		bn_mod(c[0], t0, n);
		bn_mod(c[1], t1, n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(r);
	}
}

void bn_pol_sqr_karat(bn_t *c, bn_t *a, bn_t *m, bn_t n) {
	bn_t t0, t1, a0a1;
	
		bn_null(t0);
		bn_null(t1);
		bn_null(a0a1);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(a0a1);
		
		bn_mul(a0a1, a[0], a[1]);
		bn_add(t0, a[0], a[1]);
		
		bn_pol_mul_swt(t1, a[1], m[0], n);
		//bn_mul(t1, a[1], m[0]);
		//bn_mod(t1, t1, n, u);
		
		bn_add(t1, t1, a[0]);
		bn_mul(t0, t0, t1);
		
		bn_pol_mul_swt(t1, a0a1, m[0], n);
		//bn_mul(t1, a0a1, m[0]);
		//bn_mod(t1, t1, n, u);
		
		bn_sub(t0, t0, t1);
		bn_sub(c[0], t0, a0a1);
		
		bn_lsh(c[1], a0a1, 1);
		
		bn_mod(c[0], c[0], n);
		bn_mod(c[1], c[1], n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(a0a1);
	}
}

void bn_pol_sqr_karat_old(bn_t *c, bn_t *a, bn_t *m, bn_t n) {
	bn_t t0, t1, t2;
	
		bn_null(t0);
		bn_null(t1);
		bn_null(t2);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(t2);
		
		bn_add(t1, a[1], a[0]);
		bn_sqr(t1, t1);
		bn_sqr(t0, a[0]);
		bn_sqr(t2, a[1]);
		bn_sub(t1, t1, t0);
		bn_sub(t1, t1, t2);
		
		//bn_mul(t2, t2, m[0]);
		bn_pol_mul_swt(t2, t2, m[0], n);
		bn_add(t0, t0, t2);
		
		bn_mod(c[0], t0, n);
		bn_mod(c[1], t1, n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(t2);
	}
}

void bn_pol_add(bn_t *c, bn_t *a, bn_t *b, bn_t n) {

	bn_add(c[0], a[0], b[0]);
	bn_add(c[1], a[1], b[1]);

	if (bn_cmp(c[0], n) != CMP_LT) {
		bn_sub(c[0], c[0], n);
	}
	if (bn_cmp(c[1], n) != CMP_LT) {
		bn_sub(c[1], c[1], n);
	}
}

// COUTION! Under Construction
void bn_pol_inv(bn_t *c, bn_t *a, bn_t *m, bn_t n) {
	bn_t t[2];
	bn_t r, s;
	
	bn_null(r);
	bn_null(s);
	bn_pol_null(t);
	
	TRY {
		bn_new(r);
		bn_new(s);
		bn_pol_new(t);
		
		bn_sqr(r, a[1]);
		bn_mod(r, r, n);
		bn_pol_mul_swt(r, r, m[0], n);
		bn_sqr(s, a[0]);
		bn_add(r, r, s);
		bn_pol_conj(t, a, n);
		
		//...
		
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(r);
		bn_free(s);
		bn_pol_free(t);
	}
}

void bn_pol_mul(bn_t *c, bn_t *a, bn_t *b, bn_t *m, bn_t n) {
	bn_t t0, t1, r;
	
	bn_null(t0);
	bn_null(t1);
	bn_null(r);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(r);
		
		bn_mul(r, a[1], b[0]);
		bn_mul(t1, a[0], b[1]);
		bn_add(t1, t1, r);
		bn_mul(r, a[1], b[1]);
		bn_mod(r, r, n);//
		bn_mul(r, r, m[0]);
		bn_mul(t0, a[0], b[0]);
		bn_add(t0, t0, r);
		
		bn_mod(c[0], t0, n);
		bn_mod(c[1], t1, n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(r);
	}
}


void bn_pol_mul_karat(bn_t *c, bn_t *a, bn_t *b, bn_t *m, bn_t n) {
	bn_t t0, t1, r;
	bn_t a0b0, a1b1;
		
		bn_null(t0);
		bn_null(t1);
		bn_null(r);
		bn_null(a0b0);
		bn_null(a1b1);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(r);
		bn_new(a0b0);
		bn_new(a1b1);
		
		bn_add(r, a[1], a[0]);
		bn_add(t1, b[1], b[0]);

		bn_mul(t1, t1, r);
		bn_mul(a0b0, a[0], b[0]);
		bn_mul(a1b1, a[1], b[1]);
		
		bn_mod(a1b1, a1b1, n);
		
		bn_sub(t1, t1, a1b1);
		bn_sub(t1, t1, a0b0);
		
		//bn_mul(t0, a1b1, m[0]);
		bn_pol_mul_swt(t0, a1b1, m[0], n);

		bn_add(t0, t0, a0b0);
		
		bn_mod(c[0], t0, n);
		bn_mod(c[1], t1, n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(r);
		bn_free(a0b0);
		bn_free(a1b1);
	}
}

void bn_pol_mul_karat_u_new(bn_t *c, bn_t *a, bn_t *b, bn_t *m, bn_t n, bn_t u) {
	bn_t t0, t1, r, red;
	bn_t a0b0, a1b1;
		
		bn_null(t0);
		bn_null(t1);
		bn_null(r);
		bn_null(red);
		bn_null(a0b0);
		bn_null(a1b1);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(r);
		bn_new(red);
		bn_new(a0b0);
		bn_new(a1b1);
		
		bn_set_2b(red, bn_bits(n));
		bn_mul(red, red, n);
		
		bn_add(r, a[1], a[0]);
		if (bn_cmp(r, n) != CMP_LT) {
			bn_sub(r, r, n);
		}
		
		bn_add(t1, b[1], b[0]);
		if (bn_cmp(t1, n) != CMP_LT) {
			bn_sub(t1, t1, n);
		}
		
		bn_mul(t1, t1, r);
		//bn_mod(t1, t1, n, u); //lazy out
		if (bn_cmp(t1, n) != CMP_LT) { //lazy in
			bn_sub(t1, t1, red);
		}
		
		bn_mul(a0b0, a[0], b[0]);
		bn_mod(a0b0, a0b0, n, u);
		
		bn_mul(a1b1, a[1], b[1]);
		bn_mod(a1b1, a1b1, n, u);
		
		bn_sub(t1, t1, a1b1);
		if (t1->sign == BN_NEG) {
			bn_add(t1, t1, n);
		}
		
		bn_sub(c[1], t1, a0b0);
		if (c[1]->sign == BN_NEG) {
			bn_add(c[1], c[1], n);
		}
		
		//bn_mul(t0, a1b1, m[0]);
		//bn_mod(t0, t0, n, u);
		bn_pol_mul_swt(t0, a1b1, m[0], n);
		
		bn_add(c[0], t0, a0b0);
		if (bn_cmp(c[0], n) != CMP_LT) {
			bn_sub(c[0], c[0], n);
		}
		
		
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(r);
		bn_free(red);
		bn_free(a0b0);
		bn_free(a1b1);
	}
}

void bn_pol_mul_karat_u(bn_t *c, bn_t *a, bn_t *b, bn_t *m, bn_t n, bn_t u) {
	bn_t t0, t1, r, red;
	bn_t a0b0, a1b1;
		
		bn_null(t0);
		bn_null(t1);
		bn_null(r);
		bn_null(red);
		bn_null(a0b0);
		bn_null(a1b1);
	
	TRY {
		bn_new(t0);
		bn_new(t1);
		bn_new(r);
		bn_new(red);
		bn_new(a0b0);
		bn_new(a1b1);
		
		bn_add(r, a[1], a[0]);
		if (bn_cmp(r, n) != CMP_LT) {
			bn_sub(r, r, n);
		}
		
		bn_add(t1, b[1], b[0]);
		if (bn_cmp(t1, n) != CMP_LT) {
			bn_sub(t1, t1, n);
		}
		
		bn_mul(t1, t1, r);
		bn_mod(t1, t1, n, u);
		
		bn_mul(a0b0, a[0], b[0]);
		bn_mod(a0b0, a0b0, n, u);
		
		bn_mul(a1b1, a[1], b[1]);
		bn_mod(a1b1, a1b1, n, u);
		
		bn_sub(t1, t1, a1b1);
		if (t1->sign == BN_NEG) {
			bn_add(t1, t1, n);
		}
		
		bn_sub(c[1], t1, a0b0);
		if (c[1]->sign == BN_NEG) {
			bn_add(c[1], c[1], n);
		}
		
		//bn_mul(t0, a1b1, m[0]);
		//bn_mod(t0, t0, n, u);
		bn_pol_mul_swt(t0, a1b1, m[0], n);
		
		bn_add(c[0], t0, a0b0);
		if (bn_cmp(c[0], n) != CMP_LT) {
			bn_sub(c[0], c[0], n);
		}
		
		
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t0);
		bn_free(t1);
		bn_free(r);
		bn_free(red);
		bn_free(a0b0);
		bn_free(a1b1);
	}
}


void bn_pol_mxp_basic(bn_t *c, bn_t *a, bn_t b, bn_t *m, bn_t n) {
	int i, l;
	bn_t t[2], r[2];

	bn_pol_null(t);
	bn_pol_null(r);

	TRY {
		bn_pol_new(t);
		bn_pol_new(r);

		l = bn_bits(b);
		bn_pol_copy(t, a);
		bn_pol_copy(r, t);

		for (i = l - 2; i >= 0; i--) {
			bn_pol_sqr_basic(r, r, m, n);
			if (bn_test_bit(b, i)) {
				bn_pol_mul(r, r, t, m, n);
			}
		}

		bn_pol_copy(c, r);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_pol_free(t);
		bn_pol_free(r);
	}
}


void bn_pol_mxp_slide(bn_t *c, bn_t *a, bn_t b, bn_t *m, bn_t n) {
	bn_t tab[TABLE_SIZE][2], t[2], r[2], v[2];
	bn_t u;
	int i, j, l, w = 1;
	unsigned char win[BN_BITS];

	bn_null(u);
	bn_pol_null(t);
	bn_pol_null(r);
	bn_pol_null(v);
	/* Initialize table. */
	for (i = 0; i < TABLE_SIZE; i++) {
		bn_pol_null(tab[i]);
	}
	TRY {

		/* Find window size. */
		i = bn_bits(b);
		if (i <= 21) {
			w = 2;
		} else if (i <= 32) {
			w = 3;
		} else if (i <= 128) {
			w = 4;
		} else if (i <= 256) {
			w = 5;
		} else {
			w = 6;
		}

		for (i = 1; i < (1 << w); i += 2) {
			bn_pol_new(tab[i]);
		}
		
		bn_new(u);
		bn_pol_new(t);
		bn_pol_new(r);
		bn_pol_new(v);
		
		bn_mod_pre(u, n);
		//if (bn_is_minus_one(m[0])) {
		//	bn_sub_dig(v[0], n, 1);
		//}
		//else bn_copy(v[0], m[0]);
		bn_copy(v[0], m[0]);

#if BN_MOD == MONTY
		bn_pol_set_id(r);
		bn_mod_monty_conv(r[0], r[0], n);
		bn_mod_monty_conv(t[0], a[0], n);
		bn_mod_monty_conv(t[1], a[1], n);
		//bn_mod_monty_conv(v[0], v[0], n);
		
#else /* BN_MOD == BARRT || BN_MOD == RADIX */
		bn_pol_set_id(r);
		bn_pol_copy(t, a);
#endif

		bn_pol_copy(tab[1], t);
		bn_pol_sqr_u(t, tab[1], v, n, u);

		/* Create table. */
		for (i = 1; i < 1 << (w - 1); i++) {
			bn_pol_mul_karat_u(tab[2 * i + 1], tab[2 * i - 1], t, v, n, u);
		}

		l = BN_BITS + 1;
		bn_rec_slw(win, &l, b, w);
		for (i = 0; i < l; i++) {
			if (win[i] == 0) {
				bn_pol_sqr_u(r, r, v, n, u);
			} else {
				for (j = 0; j < util_bits_dig(win[i]); j++) {
					bn_pol_sqr_u(r, r, v, n, u);
				}
				bn_pol_mul_karat_u(r, r, tab[win[i]], v, n, u);
			}
		}
		bn_pol_trim(r);
		
#if BN_MOD == MONTY
		bn_mod_monty_back(c[0], r[0], n);
		bn_mod_monty_back(c[1], r[1], n);
#else
		bn_pol_copy(c, r);
#endif
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		for (i = 1; i < (1 << w); i++) {
			bn_pol_free(tab[i]);
		}
		bn_free(u);
		bn_pol_free(t);
		bn_pol_free(r);
		bn_pol_free(v);
	}
}


void bn_pol_mxp_slide_old(bn_t *c, bn_t *a, bn_t b, bn_t *m, bn_t n) {
	bn_t tab[TABLE_SIZE][2], t[2], r[2], v[2];
	bn_t u;
	int i, j, l, w = 1;
	unsigned char win[BN_BITS];

	bn_null(u);
	bn_pol_null(t);
	bn_pol_null(r);
	bn_pol_null(v);
	/* Initialize table. */
	for (i = 0; i < TABLE_SIZE; i++) {
		bn_pol_null(tab[i]);
	}

	TRY {

		/* Find window size. */
		i = bn_bits(b);
		if (i <= 21) {
			w = 2;
		} else if (i <= 32) {
			w = 3;
		} else if (i <= 128) {
			w = 4;
		} else if (i <= 256) {
			w = 5;
		} else {
			w = 6;
		}

		for (i = 1; i < (1 << w); i += 2) {
			bn_pol_new(tab[i]);
		}
		
		bn_new(u);
		bn_pol_new(t);
		bn_pol_new(r);
		bn_pol_new(v);
		
		bn_mod_pre(u, n);
		//if (bn_is_minus_one(m[0])) {
		//	bn_sub_dig(v[0], n, 1);
		//}
		//else bn_copy(v[0], m[0]);
		bn_copy(v[0], m[0]);

#if BN_MOD == MONTY
		bn_pol_set_id(r);
		bn_mod_monty_conv(r[0], r[0], n);
		bn_mod_monty_conv(t[0], a[0], n);
		bn_mod_monty_conv(t[1], a[1], n);
		//bn_mod_monty_conv(v[0], v[0], n);
		
#else /* BN_MOD == BARRT || BN_MOD == RADIX */
		bn_pol_set_id(r);
		bn_pol_copy(t, a);
#endif

		bn_pol_copy(tab[1], t);
		bn_pol_sqr_u_old(t, tab[1], v, n, u);

		/* Create table. */
		for (i = 1; i < 1 << (w - 1); i++) {
			bn_pol_mul_karat_u(tab[2 * i + 1], tab[2 * i - 1], t, v, n, u);
		}

		l = BN_BITS + 1;
		bn_rec_slw(win, &l, b, w);
		for (i = 0; i < l; i++) {
			if (win[i] == 0) {
				bn_pol_sqr_u_old(r, r, v, n, u);
			} else {
				for (j = 0; j < util_bits_dig(win[i]); j++) {
					bn_pol_sqr_u_old(r, r, v, n, u);
				}
				bn_pol_mul_karat_u(r, r, tab[win[i]], v, n, u);
			}
		}
		bn_pol_trim(r);
		
#if BN_MOD == MONTY
		bn_mod_monty_back(c[0], r[0], n);
		bn_mod_monty_back(c[1], r[1], n);
#else
		bn_pol_copy(c, r);
#endif
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		for (i = 1; i < (1 << w); i++) {
			bn_pol_free(tab[i]);
		}
		bn_free(u);
		bn_pol_free(t);
		bn_pol_free(r);
		bn_pol_free(v);
	}
}


void bn_pol_norm(bn_t c, bn_t *a, bn_t *m, bn_t n) {
	bn_t t, r;
	
	bn_null(t);
	bn_null(r);
	
	TRY {
		bn_new(t);
		bn_new(r);
		
		bn_sqr(r, a[0]);
		bn_sqr(t, a[1]);
		bn_mod(t, t, n);
		
		//bn_mul(t, t, m[0]);
		bn_pol_mul_swt(t, t, m[0], n);
		
		bn_sub(r, r, t);
		bn_mod(c, r, n);
	}
	CATCH_ANY {
		THROW(ERR_CAUGHT);
	}
	FINALLY {
		bn_free(t);
		bn_free(r);
	}
}

void bn_pol_conj(bn_t *c, bn_t *a, bn_t n) {
	
	bn_copy(c[0], a[0]);
	if (!bn_is_zero(a[1])) {
		bn_sub(c[1], n, a[1]);
	}
	else {
		bn_copy(c[1], a[1]);
	}
}

void bn_pol_rand_in_ring(bn_t *c, bn_t n) {
	bn_rand(c[0], BN_POS, bn_bits(n)/2);
	bn_rand(c[1], BN_POS, bn_bits(n)/2);
}

void bn_pol_set_id(bn_t *a) {
	bn_set_dig(a[0], 1);
	bn_zero(a[1]);
}

void bn_pol_zero(bn_t *a) {
	bn_zero(a[0]);
	bn_zero(a[1]);
}

void bn_pol_neg(bn_t *c, bn_t *a, bn_t n) {
	if (!bn_is_zero(a[0])) {
		bn_sub(c[0], n, a[0]);
	}
	else {
		bn_copy(c[0], a[0]);
	}
	if (!bn_is_zero(a[1])) {
		bn_sub(c[1], n, a[1]);
	}
	else {
		bn_copy(c[1], a[1]);
	}
}

int bn_pol_cmp(bn_t *a, bn_t *b) {
	int c;
	
	c = bn_cmp(a[1], b[1]);
	if (c != CMP_EQ) return c;
	return bn_cmp(a[0], b[0]);
}

void bn_pol_print(bn_t *a) {

	util_print("a0: ");
	bn_print(a[0]);
	util_print("a1: ");
	bn_print(a[1]);
}

void bn_pol_reduce(bn_t *a, bn_t n) {

	bn_mod(a[0], a[0], n);
	bn_mod(a[1], a[1], n);
}

void bn_pol_trim(bn_t *a) {
	bn_trim(a[0]);
	bn_trim(a[1]);
}

int bn_pol_mul_swt(bn_t c, bn_t a, bn_t b, bn_t n) {
	
	if (b->used == 1) {
		if (bn_is_minus_one(b)) {
			bn_sub(c, n, a);
			return 1;
		}
		if (b->dp[0] == 2) {
			bn_lsh(c, a, 1);
			if (bn_cmp(c, n) != CMP_LT) {
				bn_sub(c, c, n);
			}
			return 1;
		}
		bn_mul_dig(c, a, b->dp[0]);
		if (b->dp[0] < 40) {
			while (bn_cmp(c, n) != CMP_LT) {
				bn_sub(c, c, n);
			}
			return 1;
		}
		return 0;
	}
	bn_mul(c, a, b);
	return 0;
}
